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Recent advances in realizing artificial gauge fields on optical lattices promise experimental de¬ 
tection of topologically non-trivial energy spectra. Self-similar fractal energy structures generally 
known as Hofstadter butterflies depend sensitively on the geometry of the underlying lattice, as 
well as the applied magnetic field. The recent demonstration of an adjustable lattice geometry [L. 
Tarruell et al. Nature 483, 302-305 (2012)] presents a unique opportunity to study this dependence. 
In this paper, we calculate the Hofstadter butterflies that can be obtained in such an adjustable 
lattice and find three qualitatively different regimes. We show that the existence of Dirac points at 
zero magnetic field does not imply the topological equivalence of spectra at finite field. As the real- 
space structure evolves from the checkerboard lattice to the honeycomb lattice, two square lattice 
Hofstadter butterflies merge to form a honeycomb lattice butterfly. This merging is topologically 
non-trivial, as it is accomplished by sequential closings of gaps. Ensuing Chern number transfer 
between the bands can be probed with the adjustable lattice experiments. We also calculate the 
Chern numbers of the gaps for qualitatively different spectra and discuss the evolution of topological 
properties with underlying lattice geometry. 

PACS numbers: 


I. INTRODUCTION 

The self-similar behavior of the energy spectrum of a 
charged particle moving in a lattice under a magnetic 
field has attracted great interest since its first prediction 
in 1976 [T. Such energy spectra are called Hofstadter 
butterflies due to the visually pleasing symmetry of the 
original square lattice calculation. However, experimen¬ 
tal observation of the Hofstadter butterfly requires mag¬ 
netic fields thousands of times stronger than what is cur¬ 
rently attainable for conventional solid state systems. In 
a usual solid state lattice, the length scale of the periodic 
potential is much smaller than the magnetic length. This 
obstacle has steered the efforts towards the direction of 
artificial regulation of lattice parameters or creation of 
synthetic magnetic fields. 

First attempts to experimentally realize the Hofstadter 
butterfly used superlattices which aim to enhance the lat¬ 
tice scale up to magnetic length EM]. These systems lack 
the adjustability of the lattice geometry, constituent par¬ 
ticles and interactions which are all possible with a cold 
atom system in a defect-free environment. One of the 
major accomplishments in cold atom experiments is the 
creation of effective magnetic fields through laser assisted 
tunneling M. These experiments are taking firm steps 
towards realizing degenerate Fermi and Bose gases with 
self-similar energy spectra. 

Another exciting development is the demonstration of 
an adjustable optical lattice which was used to investi¬ 
gate the dynamics of Dirac points [THE]. This lattice 
can be adiabatically changed from a dimer to a honey¬ 
comb geometry. An artificial gauge field has recently 
been added m to the same setup to realize the Haldane 
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model P2] • It is thus reasonable to expect a constant uni¬ 
form magnetic field to be simulated in this lattice. An 
adjustable optical lattice with a uniform artificial mag¬ 
netic field would make it possible to investigate not only 
the energy spectrum but also the relation between the 
lattice geometry and topological properties. 

In this paper, we study the Hofstadter butterflies and 
their topological nature for lattice geometries that can be 
created by the Zurich group. We show that this tunable 
optical lattice captures qualitatively different spectra and 
enables the study of merging of two square lattices but¬ 
terflies to form the spectrum of a honeycomb lattice. This 
merging is topologically non-trivial as an infinite number 
of gaps in the spectrum must close during this transition. 

We numerically solve the Schrodinger equation for the 
potentials reported in Ref. BE] and obtain the lowest 
two bands. In the parameter regime we are interested 
in, these bands can be accurately described by a nearest 
neighbor tight-binding approximation. We find the hop¬ 
ping strengths by fitting band structure over the whole 
Brillouin zone. Magnetic field is introduced to the tight- 
binding Hamiltonian through Peierls substitution. If the 
magnetic flux per unit cell of the lattice (in units of flux 
quantum <f>o = h/e) is a rational number >I> = p/q, en¬ 
ergy spectra can be obtained efficiently by diagonalizing 
a q x q matrix. We calculated the Hofstadter butterfly 
spectra for all the different regimes that can be accessed 
by this adjustable optical lattice. 

In the previous studies concerning adjustable or de¬ 
formed real-space lattices, the presence of Dirac points 
in the spectra has been used to mark qualitative change 
in the physics of the system. For example, a lattice that 
lacks six-fold symmetry was successfully used to mimic 
graphene El- While it is natural to expect the low en¬ 
ergy physics to be dominated by the Dirac cones if the 
Fermi surface lies close to them, it is not clear if the pres¬ 
ence of Dirac cones qualitatively determine the global 
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shape of the energy spectrum. In this sense, one can ask 
whether the presence of Dirac cones at zero magnetic 
field for two different systems necessitates their topologi¬ 
cal equivalence at non-zero magnetic held. We show that 
this equivalence is not guaranteed, be. two lattices both 
of which have Dirac points at zero held may have quali¬ 
tatively different behavior at non-zero held and may have 
different topological characters. More specifically, as the 
lattice alters from the checkerboard to the honeycomb 
geometry, the energy spectrum changes from two sepa¬ 
rate square lattice butterflies to the honeycomb lattice 
butterfly. However, the energy gap around zero closes at 
different points during the transition depending on the 
denominator q of the hux $ = p/q. We characterize the 
topological nature of the transition by calculating the 
Chern numbers for the major gaps of the spectrum. 

The paper is organized as follows: In the next section, 
we summarize the properties of the lattice generated by 
the Zurich group and outline our numerical procedure for 
the tight-binding description. The third section, we dis¬ 
cuss the effect of the magnetic field and solve the result¬ 
ing Hamiltonian to obtain the Hofstadter butterfly in the 
interesting parameter regime. Sec. IV presents the cal¬ 
culation of the first Chern numbers in the butterfly and 
a detailed discussion of the transition. We conclude in 
Sec. V with a summary and experimental consequences 
of our results. 


II. THE MODEL 

We start by considering the single-particle spectrum 
in a two-dimensional periodic potential in the absence 
of the magnetic field. As long as the Wannier functions 
for the lowest two bands are well localized, the effect 
of the magnetic field can be included by modifying the 
hopping parameters via the Peierls substitution |l3j . To 
obtain the tight-binding parameters, we first calculate 
the energy spectrum of the system by numerical solution 
of the Schrodinger equation. The real-space Hamiltonian 
without the magnetic held is 

K 2 ( d 2 d 2 \ 

w = -^(ae + v) +F(l - !,) ' (1) 

where we use the potential V(x,y) reported by the Zurich 
group [7], 

Q 

V(x,y) = — 14 cos 2 (kx + -) — V x cos 2 {kx) — V y cos 2 {ky) 

— 2a\JV x V y cos (kx) cos (ky) cos(^), (2) 

with k = 27t/A, A being the wavelength of laser beams 
creating the optical lattice, a is the visibility of the in¬ 
terference pattern. The phases 6 = n and tp = 0 are set 
to the experimental values. The potential depths we use 
are in the experimentally reported regime. The confine¬ 
ment along the ^-direction is strong so that the system 
consists of independent layers of two-dimensional sheets. 


As the parameters 14,14 and V y are varied, the Bra- 
vais lattice defined by the periodicity of the potential 
does not change. However, within a unit cell, the num¬ 
ber and positions of potential minima vary determining 

the overall lattice geometry. If the potentials are deep 

, 2 

enough compared to recoil energy E R = " A2 , states 

localized around these minima dominate the low energy 
physics and determine the structure of the lowest bands 
in the spectrum. The lattice formed by these states span 
a range of geometries varying from the checkerboard to 
the honeycomb lattice. Within this range, the most im¬ 
portant qualitative change is the creation of Dirac points 
in fc-space, which was investigated in detail both exper¬ 
imentally and theoretically mmm- At zero magnetic 
field, we follow a similar procedure to obtain a tight- 
binding description. 

We rotate the coordinate system by 7t/4 and then use 
X/V2 to make positions x,y dimensionless making the 
unit cell a one-by-one square in the x — y plane. All the 
potentials and energies are measured in units of recoil 
energy E R . We obtain numerical solutions of this dimen¬ 
sionless Hamiltonian for the lattice parameters defined 
above. We discretize the real-space Hamiltonian Eq[l] 
and numerically diagonalize the resulting energy matrix 
for 2500 points per unit cell with boundary conditions 
determined by the Bloch theorem. 

For the parameter regime we are interested in, the lat¬ 
tice potentials are chosen deep enough to create nar¬ 
row and symmetric energy bands, as the correspond¬ 
ing Wannier functions are well localized in real-space. 
Thus, a bipartite symmetric tight-binding approxima¬ 
tion to these bands successfully describes the system. 
We construct the lattice with Bravais lattice vectors 
Rm 1 ,m 2 = JTiidi +771202 for (mi,m 2 ) £ Z for dimension¬ 
less a[ = x, a ~2 = y. There can be one or two potential 
minima within the unit cell depending on the parameters. 
As we are interested in the emergence of Dirac points 
which corresponds to the transition from the dimer to 
the honeycomb-like lattice, we choose a unit cell with a 
two-state basis (see FigjI]). These states, labeled as A 
and B, are localized at 

dA. = ^(~x + y), d B = ^=(x-y), (3) 

where d is the distance between the sites A and B in units 
of the lattice constant. In the experiments, adjusting the 
potential parameters modifies this distance d which ef¬ 
fectively changes the hopping parameter between sites 
A and B. That is to say, without the magnetic field, 
the distance d is embedded into tunneling amplitudes 
to,ti in the tight-binding Hamiltonian. As we work in 
the regime where the lowest two bands are narrow, we 
limit out tight-binding description to nearest-neighbor 
hopping. Each site has three nearest-neighbors; one in 
the same unit cell and two outside. Each state of type A 
is connected to only states of type B. The hopping pa¬ 
rameter t 0 connects the states within the unit cell while 
t\ links neighboring unit cells. Due to the bipartite sym- 
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FIG. 1. (Color online) The real-space potential and corre¬ 
sponding tight-binding lattice. The unit cell remains square 
for all experimental parameters [8]. The tunneling amplitudes 
to,ti are calculated by a fit to numerically obtained two low¬ 
est band energies. For V x = 14,14= 0.79, V y = 6.45 in units 
of Er and a = 0.9, we find to = 0.0391, t\ — 0.0361. Within 
a unit cell, two well localized Wannier functions at poten¬ 
tial minima (A,B) form the basis for tight-binding treatment. 
Modification of the potential depth parameters changes the 
distance d between the sites A and B altering the tunneling 
rates between these sites. 

metry of the lattice, the lowest two bands are symmetric 
in energy. The tight-binding Hamiltonian follows: 

R = — ^ m 2 , A)(ml, m2, B\ 

mi ,7722 

+ ii|mi + 1, m 2 , A)(ml, m2, B | 

+ f 1 17T?.i, m 2 — 1, A)(ml, m2, B\ + h.c. j, 

( 4 ) 

where h.c. stands for Hermitian conjugate. The self en¬ 
ergy term eo is excluded so that the energies are symmet¬ 
ric around zero. 

The tight-binding Hamiltonian is readily diagonalized 
by Fourier transformation 

e (toRi) = 

± y4jj + 2 t\ + 2f 0 ii(cos k x + cos k y ) + 2t\cos(k x — k y ), 

( 5 ) 


eluded) reduce to just two parameters in the tight- 
binding model to, U ■ When the energy is scaled by a 
tunneling rate, just a single parameter ti/to determines 
the spectrum. Each different lattice is represented by 
its t\/to ratio which gives a different butterfly spectrum 
when the magnetic field is applied. 

As we only include the nearest-neighbor hoppings, the 
lattice is bipartite. In order to have a Dirac point, 
Schrodinger equation must have a zero eigenvalue. For a 
state with zero energy the total incoming current to any 
lattice site must be zero. Each site has three neighbors 
all of which reside on the same sublattice. Thus, the wave 
function on each neighbor can vary by at most a phase 
factor. The existence of a zero energy solution depends 
on the solution of to + + e 1 ^ 2 ) = 0. Therefore, we 

obtain the condition for the existence of Dirac points as 
R > 4 at zero magnetic field in agreement with Ref. Il4j. 
We stress that the existence of Dirac point at zero mag¬ 
netic field does not imply the existence of zero energy 
solutions at finite field. 

The effect of an external magnetic field can be included 
in the continuum Hamiltonian [T] by minimal coupling 
p —> p — eA. However, the modification of the lowest 
few bands by the magnetic field can be accurately de¬ 
scribed through tight-binding with complex hopping pa¬ 
rameters. This procedure, called the Peierls substitution, 
was found to be precise as long as the tight-binding de¬ 
scription at zero held holds [15] . For the recent cold atom 
experiments, this substitution method also captures the 
physical process mm- The artificial magnetic fields in 
optical lattices are simulated by modification of the hop¬ 
ping between neighboring sites. 

We use the Landau gauge A = (0, B 0 x, 0) for the mag¬ 
netic held B = B 0 z. The tunneling amplitudes acquire a 
complex phase 

t l m,n\n,j)(m,i\ -> e- lQ ™t l ^ n \n,j)(m,i\, (6) 

where | f^i) ^ ^ f° r l a tti ce points m = 

(mi, m 2 ), n = ( 711 , 712 ) and sites (i, j). The hopping phase 
between two arbitrary lattice points is: 


where k x (k y ) is the quasi momentum along x(y)~ 
direction. By htting this expression to the numerical 
results (Eq§, we find the tunneling amplitudes to,ti 
together with the self energy eo for any optical lattice 
potential in an efficient way. In the earlier studies con¬ 
cerning the merging of Dirac points, the hopping parame¬ 
ters were determined by htting these two energies specif¬ 
ically around the Dirac cones nu As we are interested 
in the global properties of the spectrum, we minimize 
the least-square error in energies of both bands over the 
entire Brillouin zone. For our calculations, this error is 
limited by 4.5% as an average with respect to the band 
widths and reduces as we move towards the honeycomb 
geometry. 

Therefore, the three parameters from the experiments 
(six parameters if the phases and the visibility are in- 


Q ij = 2tt$ 

m n ^ ' 1 


' (Rn Rm 4 “ dj di) 

Rn + Ri 


X ■ 


dj + di ) 


Therefore, the tight-binding Hamiltonian under the mag¬ 
netic held becomes 


H = - {t 0 \mi,m 2 ,A)(mi,m 2 ,B\e l2n ® dm i 

mi ,rri2 

+ ti |mi + 1, m 2 , A)(mi, m 2 , B\e~ lM ^ + ^ 

+ ti\mi,m 2 — 1, A)(m\, m 2 , B\e~ i27T<s?( - d ^ 1 ) T ™ 1 -f h.c. 

( 7 ) 


where $ = (j)/4>o stands for the total magnetic hux <f> 
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passing through each unit cell over flux quantum (f >o = 
h/e and d' = + . 

In the next section, we obtain the eigenvalues of this 
Hamiltonian and analyze the transition from the checker¬ 
board to the honeycomb lattice. 

III. THE TIGHT-BINDING SOLUTION 


Since we employ the Landau gauge with the vector po¬ 
tential A along one of the primitive vectors of the lattice, 
we retain the translational symmetry along one direction. 
Thus, the m 2 dependent part of the wave function can 
be separated 

= e ikyTn2 e ~ in ' s ’ d ' mi f mi , (12) 


To find the eigenvalues, we consider rational flux values 
$ = p/q where p and q are co-prime integers. Eigenstates 
of the Hamiltonian EqjT] can be expanded in real-space 

+> = 55 ^t u rn 2 \ m ^ m 2,A)+^ u m 2 \rni,rn 2 ,B). 

777 1,777-2 

( 8 ) 

The Schrodinger equation for the coefficients at sites A 
and B reads 


+mi,m 2 = - h eM{mi+1/2) Tpt+l,r 

c i 2 n®d,'m,! n i A 

U m\ ,777-2 


_ —i 2 n$(l—d')mi i A 

r mi ,m 2 —1 ’ 


(9) 


fl lA — _ f p -i2n$d' (ttii-1/2) /B 

^7711,7712 _ C l e 1 777-1 — 1,71 




$ 


B 

777-1 ,7772 


^ g^27T^«J> (1 — d! ) 7771 


■0 


B 

777-1,7772 + 1 ’ 


( 10 ) 


where e is the eigenvalue. These two sets can be combined 
to give a quadratic equation 

+lm 2 = (io + 2 ^)V’m 1 ,m 2 

+ t 0 ti{e i “ 1 <,m 2 -i+e-^ d '<_ 1 ,m 2 

, Slim! ii , i7r$d' / A \ 

' ^ +7771,7772 + 1 ^ ^ +777l+1,7772 J 

I J.2 f —i27T$(7771 -1) -77r$d' /A 

' 6 1 1 ^ C +7771-1,7772-1 

+ e i “ 1 e i ™'< +1 , m2+1 }. (11) 


where k y is the quasi momentum along the y-direction. 
The second exponent term is introduced to make the 
eigenvalue equation independent of d' through a gauge 
transformation. Reorganizing Eq[TT] with the help of 
Eq[l2] results in a one-dimensional difference equation, 


A f mi = (1 + -e- i(2,r$ ( mi - 1)+fc ^ 

V *o 

+ 2 cos(27r$mi + k y )f mi 

+ (l + / mi+1 

with eigenvalues 



(13) 


A = 


(*g + 2tj) 

toti 


(14) 


This equation is periodic with q and the Bloch condi¬ 
tion then allows us to write 


fm 1+q {k x ) = e i ^fm 1 {kx). (15) 

The Bloch condition reduces the calculation of the eigen¬ 
values to diagonalization of a q x q matrix at every value 
of k x ,k y . Because of the 9 -fold periodicity, unique val¬ 
ues of k x are limited to the region — ^ < k x < ^ while 
—7 r < k y < 7 r defining the magnetic Brillouin zone. Ev¬ 
ery value of A found at k x ,k y defines two e symmetric 
around zero reflecting the bipartite symmetry. Here, A 
are the eigenvalues of 


2cos ( 271 -10 + 


1 -|_ A e i ( 27r f o+ky) 
to 


l 1 A e -i ( 2 ”'? 0+k y 
to 


e i 9 fc x _|_ *i e -d 27r q- k y -qkA 


2cos(27r|l + fc y ) 1 + |e i( 27 r f 1+ ^) 


e -i 9 fe x _|_ *i e i(27r|-fc B --3fc x )A 


1 ti e i(2ir^(q-2)+k y ) 

to 


to 


1 + k e -i( 2 -f ( ?- 2) +M 2cOs(27rj|(g-l) + fc„) 

(16) 


While it is possible to do numerical diagonalization at each point in the Brillouin zone, such a numerically costly 
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approach is not necessary as long as one is interested in 
the width of the bands and gaps in the spectrum. The 
maxima and the minima of the bands appear at the same 
k -space points for all bands at a given flux. Determina¬ 
tion of these special points enables us to plot the spec¬ 
trum at any flux by diagonalizing two q x q matrices 

[UHSIII7]. 

The secular determinant of the eigenvalue equation has 
a simple dependence on k x . k y . Due to the 27 r /<7 period¬ 
icity of the energy in both k x and k y , there is only one 
momentum-dependent term, H(k x ,k y ), in the determi¬ 
nant. This term is independent of A and can be calcu¬ 
lated as 

D (k x , k y ) = (- 1) 9+1 cos (qk x ) + (-1 cos (qk y ) 

+ (-l ) p(<? - 1)+9+1 V cos (q(k x + ky)). 

(17) 

The extrema of this expression define the points in the 
Brillouin zone at which all bands have their maximum or 
minimum energy. These critical k points can be found as 
given in Table |T| 

Both the calculation of the tight-binding parameters 
for the experimental optical lattice and the calculation 
of the Hofstadter butterfly for those tunneling rates are 
numerically efficient. This allows us to calculate the but¬ 
terfly spectra in all the experimentally accessible regimes. 
In the next section, we focus on qualitatively different 
spectra that can be obtained and discuss the non-trivial 
transition between them. 


IV. THE ENERGY SPECTRA 

The experiment in Ref. [7] demonstrates an extremely 
versatile optical lattice. While it is possible to calculate 
the butterfly spectrum for any set of experimental pa¬ 
rameters, in this paper, we are mainly concerned with 
what kinds of topologically distinct spectra are accessi¬ 
ble. The Hofstadter butterfly features an infinite number 
of gaps due to its fractal nature and closure of any one 
of these gaps creates a topologically different spectrum 
in principle. However, the smaller gaps in the spectrum 
are unlikely to be observed experimentally due to the fi¬ 
nite size and temperature effects [TSj. The closing of the 
largest gaps in the spectra signals drastic changes in the 
behavior of the system and is more likely to have exper¬ 
imental signatures. 

As characterized by their major gaps, the butterfly 
spectra for the adjustable lattice generally falls into one 
of the two categories. The spectrum is either similar to 
two square lattice Hofstadter butterflies separated by a 
large gap or it is similar to the butterfly of the honey¬ 
comb lattice (Hofstadter-Rammal butterfly). As we show 
below, the transition between these two general classes 
is not abrupt. The major gap near zero energy does not 
close at a single critical lattice parameter, but, this gap 


fractures into smaller gaps through gap-closings at prin¬ 
ciple fractions of flux. We discuss the physics of this 
process by identifying four representative points in the 
parameter space and displaying their spectra m- 

The geometry of the adjustable lattice in Ref. [7] 
changes mainly by the reorganization of the potential in¬ 
side a unit cell while the Bravais lattice always remains 
square. Hence, the different regimes such as ’’dimer”, 
’’checkerboard”, etc. can be best understood by investi¬ 
gating the potential within a real-space unit cell. This 
potential can have one or two minimum points and the 
transition from one to two minima signals the change 
from the checkerboard to the dimer lattice. If the po¬ 
tentials are deep enough to form a tight-binding lattice, 
the Wannier functions localized around these minima 
have very small overlaps with the neighboring unit cell. 
The lowest two states form rather narrow bands which 
are separated by a gap controlled by the overlaps inside 
the unit cell. Thus, within the dimer and checkerboard 
regimes, the spectrum under the magnetic field consists 
of two square lattice Hofstadter butterflies. 

As an example, in Figj2] we display the results of 
our calculations for the experimental parameters V x = 
4.5,14 = 1 ,V y = 5 and a = 0.95. For these param¬ 
eters, the real-space unit cell has two minima close to 
each other but far away from the minima in the neighbor¬ 
ing cells. This arrangement was called the dimer lattice 
(see Fig.l in Ref.[?]). We find the tunneling amplitudes 
t 0 = 0.3488, t\ = 0.0843 for these parameters and calcu¬ 
late the corresponding butterfly. The spectrum is clearly 
formed by two square lattice Hofstadter butterflies sepa¬ 
rated by a large gap. This separation is controlled by the 
in-cell hopping parameter to- The upper (lower) butterfly 
is formed by the square lattice of anti-bonding (bonding) 
states in each unit cell. 

The lattice gets closer to the honeycomb regime as 
the two minima within the unit cell move away from 
each other. The overlap between localized states at each 
minimum decreases with the increasing distance between 
them. Simultaneously, the distance between minima 
in adjacent cells decreases, increasing the corresponding 
tunneling rate. While the tight-binding amplitudes are 
not determined solely by the distance between the min¬ 
ima, this rough picture points to to decreasing and t\ 
increasing as honeycomb regime is approached. The zero 
field band structure shows an abrupt transition evidenced 
by the creation of Dirac points at the zone corners. As 
calculated above, the Dirac points exist for | and 
this value was generally used to mark the boundary of 
the honeycomb regime 13 El. 

We calculate the spectrum at non-zero field for lattices 
which have Dirac points at zero field. We find that the 
existence of the Dirac points at zero field does not guar¬ 
antee the existence of zero energy states at finite mag¬ 
netic held. As a concrete example, in Fig(3]we plot the 
butterfly for t 0 = 0.0691, t\ = 0.0361 which corresponds 
to experimental values V x = 12, V x = 0.79, V y = 6.45 and 
a = 0.9. For these values, the major gaps of the spectrum 
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TABLE I. Critical values for extremum points of the energy bands for flux = |. 
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FIG. 2. (Color online) Energy spectrum as a function of flux 
per unit cell. For 14 = 4.5, V x — 1, V y = 5 and a = 0.95, 
creating a lattice in the Dimer regime. All energies are in 
units of Er, flux is measured in units of flux quantum. Tun¬ 
neling amplitudes are to = 0.3488 and t\ = 0.0843 in units of 
recoil energy. The spectrum is formed by two square lattice 
Hofstadter butterflies separated by a large gap controlled by 
to- The width of both butterflies are proportional to ti. In¬ 
set is used to display energy bands as a function of k x , k y at 
zero magnetic held. The gap at zero magnetic held is mostly 
retained at non-zero held. 


resemble two square lattice butterflies touching only at 
their edges. The gap between these two segments lies un¬ 
broken for all flux values. Hence, although the zero field 
spectrum has Dirac points like the honeycomb lattice, 
the global shape of the spectrum under magnetic field is 
closer to two stacked square lattice butterflies. 

As another example, consider a lattice deeper in the 
honeycomb regime. Tunneling amplitudes are to = 
0.0441, ti = 0.0361 which can be realized by 14 = 
13.5,14 = 0.79, V y = 6.45 and a = 0.9. Again, although 
Dirac points exist at zero field, there are large gaps at 
zero energy. However, the overall central gap observed in 
the Fig|3]is now broken up into four parts. These closings 
happen at flux p/q = 1/2,1/3,2/3 but nowhere else. 

The critical ratio ti/to for which the central gap closes 
at a particular flux can be calculated analytically. If the 
gap is closed, the Eq|9) Eq[l0] must admit a zero energy 


FIG. 3. (Color online) Hofstadter butterfly for 14 = 12,14 = 
0.79,14 = 6.45 and a = 0.9, which is at the boundary of the 
honeycomb regime. Tunneling amplitudes are to = 0.0691 and 
ti = 0.0361 in units of recoil energy. Comparing with FigH 
one can see that the two stacked Hofstadter butterflies get 
closer and touch only at the edges. As displayed in the inset, 
there are Dirac points at zero magnetic field band structure. 
However, for all non-zero magnetic fields, there is a finite gap 
at zero energy. 


solution. At zero energy, these equations relate the coef¬ 
ficients of the wave function in adjacent unit cells, q fold 
iteration connects two coefficients which are equivalent 
up to a phase due to the Bloch condition: 


(_1)p0? 


-D P ^ q = ( _ 1} 


9g ik y<l 


(18) 


This relation is satisfied for some point in the Brillouin 
zone iff 


h>^_ 

t 0 ~ 2 1 /?' 


(19) 


Thus, for any flux 4> = |, there is no gap near zero energy 

if 4 > regardless of the value of p. As the lattice 
evolves from the dimer to the honeycomb geometry, the 
first point where the central gap closes is at zero magnetic 
field for t\ = which is the formation of Dirac points 


at zone corners. However. 


for 0.5 < ± < -4 ~ 0.7071, 

£ 0 v 2 
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FIG. 4. (Color online) Hofstadter butterfly for 14 = 
13.5,14 = 0.79, V y = 6.45 and a = 0.9, deeper in the 
honeycomb regime. Tunneling amplitudes are to = 0.0441 
and t\ = 0.0361 in units of recoil energy Er. The central 
gap is closed at the edges (4> = 0,1) and the flux values 
p/q = 1/2,1/3, 2/3 but nowhere else. For ^ = 0.8186, Eq. 


19 


shows that all fluxes with q > 4 have a gap at zero energy. 
The central gap is now split into four parts. Comparing the 
inset with the almost identical inset of Fig[3] reveals the in¬ 
ability of zero magnetic field band structure to determine the 
spectrum for non-zero field. 


the gap at zero energy remains open at any non-zero 
magnetic field. The next flux for which the gap closes 
is < I ) = 2 , followed by $ = | at = ^73 ~ 0.7937. 

As an example, we plot the butterfly for = 0.8186 in 
Fig|4] where one can observe the existence of the central 
gap except for three points discussed above. 

The stitching of two square lattice butterflies to form a 
honeycomb lattice butterfly is a topologically non-trivial 
process. As we calculated above, this evolution involves 
the creation and then subsequent closing of infinitely 
many gaps. Even more surprising is the fact that the 
critical values for gap-closings at zero energy follow a 
simple yet non-analytic function which depends only on 
the denominator of the flux. The vanishing of the central 
gap is completed only when t\ = to, which identically 
corresponds to the honeycomb lattice. It is thus natu¬ 
ral to ask whether the lattice geometry that was used 
to simulate graphene [Hj retains this similarity at finite 
field. In Fig(5j we show the butterfly for the experimental 
parameters used to simulate graphene. For these param¬ 
eters, we calculate that the central gap is closed for all 
flux values with denominators less than q = 9. However, 
comparing this spectrum with the butterfly for an ideal 
honeycomb lattice FigjGj we see that while the gaps at 
zero energy have not completely closed, they are suffi¬ 
ciently small so that the overall shape of the butterfly is 
almost the same. As long as the energy resolution of a 
measurement is coarser than the size of these gaps, the 
same parameters can be used to simulate graphene under 


FIG. 5. (Color online) Hofstadter butterfly for 14 = 14, 14 = 
0.79, V y = 6.45 and a = 0.9, the parameters used to simulate 
graphene 0. Tunneling amplitudes are to = 0.0391 and t\ = 
0.0361 in units of recoil energy. For these values, the spectrum 
nearly matches the ideal honeycomb butterfly Fig[6j except 
for minor gaps. The central gap is open for the flux values 
with denominator q > 9, producing narrow gaps near zero 
energy. As seen in the inset, the change in the zero field band 
structure is only due to the movement of Dirac points. 



FIG. 6. Hofstadter butterfly of the ideal honeycomb lattice 
(to = ti) as obtained in Ref. [T 6 . Energy is plotted in units 
of to- Notice that there is no gap at zero energy for any 
value of flux in contrast to Figj5] Still, the overall structure 
is remarkably similar. 


a magnetic field. 

While the rough structure of the spectrum as charac¬ 
terized by the major gaps follows the physical picture 
given in this section, it is important to describe topo¬ 
logical properties of the evolution in detail. In the next 
section, we calculate the Chern numbers for the gaps in 
the spectrum and show that bands exchange topologi- 










































cal invariants through closing and subsequent opening of 
gaps. While infinitely many closings take place, all of 
them occur for minor gaps in the spectrum and are un¬ 
likely to be observed within current experimental energy 
resolution. 


V. CHERN NUMBERS 

The energy bands for a periodic potential under a mag¬ 
netic field have a topological invariant called the first 
Chern number which is directly related to the contribu¬ 
tion of that band to Hall conductivity m- As long as 
a band is isolated by gaps in energy, this integer can¬ 
not change due to changes in the lattice potential or the 
magnetic field. Thus, integer Chern numbers can be used 
to label the gaps in the butterfly spectra [2D- If the 
Fermi energy lies in a gap labeled with Chern number 
i/, the Hall conductivity for a non-interacting system of 
fermions is a xy = ve 2 /h. 

Since Chern numbers are topologically protected 
against small perturbations and directly related to mea¬ 
surable quantities, they have attracted both theoretical 
and experimental interest recently [22HMI- A number of 
methods for Chern number and Berry curvature measure¬ 
ment in cold atom systems have been proposed and im¬ 
plemented. A method which is not only useful for exper¬ 
imental Chern number measurement but also for numer¬ 
ical calculation relies on the Streda formula G3- Streda 
formula connects the response of the density to changes 
in magnetic flux with the Hall conductivity, hence the 
Chern number. 

_ dN 
° xv ~ dB' 

where N is the number of levels lying below the Fermi 
energy. As density probes are more common in cold atom 
experiments compared to transport measurements, this 
relation is exceptionally suited for ultracold systems [Si- 
Adapting this formula to the system in consideration, the 
Chern number for any gap is 


where n is the density per unit cell for fermions with 
Fermi energy lying in that gap. Numerically we calculate 
Chern numbers by implementing the derivative by a finite 
difference. The calculation is reduced to counting of the 
eigenvalues of the matrix Eq |16| below the Fermi energy 
for two closely spaced flux values [IS . 

We first test the observation made in the previous 
section that the largest gaps except the central gap in 
the spectrum remain open during the evolution from the 
checkerboard to the honeycomb lattice. As can be seen 
in Fig{7] calculated for flux $ = 1/3, the gaps with Chern 
numbers ±1 remain open for all 0 < ~ < 1. In one limit, 
these gaps can be interpreted as the major gaps of two 
stacked square lattice butterflies while in the other limit 
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FIG. 7. (Color online) The evolution of bands from the 
checkerboard to the honeycomb regime as a function of ^ 
for flux <1? = 1/3. Gaps are labeled with their Chern num¬ 
ber. The gaps with non-zero Chern numbers are open for 
all 0 < < 1. The central gap closes at = ^ 373 - Near 

= 0 . 1 , the gaps are the gaps of the two stacked square lat¬ 
tice butterflies. In the other limit, they connect to the largest 
gaps of the honeycomb lattice without gap closings. 


they are the largest gaps of the honeycomb lattice. The 
central gap with Chern number zero closes at |^ = ■ 

The simple evolution for the largest gaps is not re¬ 
peated by smaller gaps. In general, many small gaps are 
closed and reopened as increases m- The number of 
bands at flux $ = p/q is 2 q for two stacked square lattices 
and 2q— 1 for the honeycomb lattice. Apart from the clos¬ 
ing of the central gap, the number of gap closures must 
be equal to the number of gap openings. However, two 
bands can change their Chern numbers through touching 
and splitting, conserving the total Chern number. 

As an example, we display the bands at $ = I as a 
function of in Fig[ 8 l Near cs 0 , the gaps in the 
spectrum have the Cnern numbers 1, 2, -2, -1, 0, 1, 2, 
-2, -1 in the order of increasing energy. This sequence is 
replaced by 1, 2, 3, -1, 0, 1, -3, -2, -1 in the honeycomb 
lattice limit I 1 = 1. We see that in addition to the central 
gap closing at = ^ 73 -, the third (eighth) and the fourth 
(seventh) bands touch and exchange Chern numbers at 
/l = 0.5. 

In general, for larger denominator q, there will be 
more contact between the bands throughout the evolu¬ 
tion Fig[9| As the number of bands for the honeycomb 
lattice is less than the number of bands for two stacked 
square lattice butterflies by just one, all the gaps except 
the central gap are related by a finite number of Chern 
number transfers between the bands. In that sense, it is 
not surprising that the Diophantine equation which con¬ 
straints the gap Chern number is the same for square and 
honeycomb lattices ESI . 

In particular, the experimental setup used in Ref. [5] to 
mimic graphene will differ from graphene under a mag- 
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FIG. 8 . (Color online) The evolution of bands from the 
checkerboard to the honeycomb at $ = 1/5 as a function of 
. Gaps are labeled with their Chern number. The sequence 
of the Chern numbers in the spectrum near ^ = 0.1 is 1, 2, 
- 2 , - 1 , 0 , 1 , 2 , - 2 , -1 as expected for two stacked square lattice 
butterflies. In the other limit, the sequence is 1, 2, 3, -1, 0, 1, 
-3, -2, -1 valid for the honeycomb lattice. The change in the 
sequence is possible due to the third (eighth) and the fourth 
(seventh) bands touching and exchanging Chern numbers at 
= 0.5. The central gap closes at ^ = ^375 ■ 


netic field for flux values with denominators greater than 
q = 8. The largest gap that is open for the optical lat¬ 
tice but not for the ideal honeycomb lattice has an energy 
width which is almost one percent of the full band width. 
Unless the experiment probes with this energy resolu¬ 
tion, the butterflies of the optical lattice and graphene 
are equivalent. 

As the experiment can access all A values, it is possible 
to observe adiabatic Chern number transfer between the 
bands. Consequently, a lattice with large enough gaps 
which are adiabatically connected to honeycomb lattice 
gaps can be manufactured. A measurement of the inte¬ 
ger Chern numbers within the largest gap sequence near 
zero energy is possible with the experimental parameters 
used in Ref. [5]. Such a measurement would amount to 
a demonstration of the integer quantum Hall effect for 
Dirac fermions |29j in cold atom systems. 


VI. CONCLUSION 

Realization of artificial magnetic fields in optical lat¬ 
tices enables the observation of Hofstadter butterfly spec¬ 
trum with cold atoms. This self-similar energy spectrum 
is highly specific to the lattice geometry. Hence, the 
latest advances in creating adjustable lattice geometries 
provide the opportunity to study the dependance of the 
Hofstadter butterfly to the lattice parameters. 

We calculated the energy spectra for experimentally 
realized potentials [7l at non-zero magnetic field. The 
effect of the magnetic field is taken into account by using 



FIG. 9. (Color online) Evolution of the bands and gap Chern 
number for = 1/13, $ = 5/13. For large denominators, the 
evolution contains many band crossings and Chern number 
transfers. These transfers are necessary for the sequence of 
gap Chern numbers to connect the two limits. Notice that 
the largest gaps remain open throughout. Chern numbers for 
the smallest gaps are not displayed for visual clarity. 
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Peierls substitution on a nearest neighbor tight-binding 
model for the potential. 

Our results indicate that if an artificial magnetic field 
is created in the tunable lattice used in Ref. it would 
be possible to observe the evolution of the spectrum from 
two square lattice butterflies to a honeycomb lattice but¬ 
terfly. As the largest gaps remain open throughout the 
transition, Chern numbers measured at any point within 
these gaps are equivalent to the ideal honeycomb lattice 
spectrum. The most prominent gap that closes during 
the transition is the central gap, which vanishes by se¬ 
quential creation of Dirac points in the spectrum at ra¬ 
tional flux values. We calculated these critical points for 
any rational flux $ = | in terms of the tight-binding 

parameters ^^-. 

Our most important conclusion is that the existence 
of Dirac points at zero magnetic field does not guarantee 
that the Hofstadter butterfly is similar to the spectrum 
of the ideal honeycomb lattice. As the A ratio is in¬ 
creased, the general structure of the butterfly spectrum 
approximates the honeycomb butterfly better. The two 
butterflies are identical only when A = 1. For any ratio 


smaller than one, infinitely many but extremely narrow 
gaps are present. Bands in the spectrum exchange Chern 
numbers through the closing and reopening of these mi¬ 
nor gaps. 

In principle, the evolution of this lattice makes it pos¬ 
sible to observe Chern number transfer between the ad¬ 
jacent bands through band touching. However, we found 
that the gaps which close and reopen are quite narrow 
compared to the full band width. The energy resolution 
to investigate this process can only be obtained by using 
large lattices so that finite size effects do not wash out 
these small gaps. If these small gaps are not resolved, the 
artificial graphene created by the Zurich group would suc¬ 
cessfully simulate the Hofstadter butterfly of graphene. 
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